arXiv:1507.05063v2 [math.NA] 31 Aug 2016 


Noname manuscript No. 

(will be inserted by the editor) 


An extrapolation cascadic multigrid method combined with a 
fourth-order compact scheme for 3D Poisson equation 


Kejia Pan • Dongdong He • Hongling Hu 


Received: date / Accepted: date 


Abstract Extrapolation cascadic multigrid (EXCMG) method is an efficient multigrid method 
which has mainly been used for solving the two-dimensional elliptic boundary value prob¬ 
lems with linear finite element discretization in the existing literature. In this paper, we 
develop an EXCMG method to solve the three-dimensional Poisson equation on rectangu¬ 
lar domains by using the compact finite difference (FD) method with unequal meshsizes in 
different coordinate directions. The resulting linear system from compact FD discretization 
is solved by the conjugate gradient (CG) method with a relative residual stopping criterion. 
By combining the Richardson extrapolation and tri-quartic Lagrange interpolation for the 
numerical solutions from two-level of grids (current and previous grids), we are able to 
produce an extremely accurate approximation of the actual numerical solution on the next 
finer grid, which can greatly reduce the number of relaxation sweeps needed. Additionally, 
a simple method based on the midpoint extrapolation formula is used for the fourth-order 
FD solutions on two-level of grids to achieve sixth-order accuracy on the entire fine grid 
cheaply and directly. The gradient of the numerical solution can also be easily obtained 
through solving a series of tridiagonal linear systems resulting from the fourth-order com¬ 
pact FD discretizations. Numerical results show that our EXCMG method is much more 
efficient than the classical V-cycle and W-cycle multigrid methods. Moreover, only few CG 
iterations are required on the finest grid to achieve full fourth-order accuracy in both the 
L 2 -norm and L“-norm for the solution and its gradient when the exact solution belongs to 
C 6 . Finally, numerical result shows that our EXCMG method is still effective when the ex¬ 
act solution has a lower regularity, which widens the scope of applicability of our EXCMG 
method. 
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1 Introduction 

Poisson equation is a partial differential equation of elliptic type with broad application in 
electrostatics, mechanical engineering, theoretical physics and geophysics. The Dirichlet 
boundary value problem for the three-dimensional (3D) Poisson equation has the following 
form: 

f u xx + Uyy + u 7Z = f{x, y, z), in Q, 
u(x,y,z)= g(x, y,z), ondO, 

where Q is a 3D rectangle domain and dQ is its boundary. Here we assume that the forcing 
function f(x, y, z), the boundary function g(x, y, z) and the exact solution u(x, y, z) are con¬ 
tinuously differentiable and have the necessary continuous partial derivatives up to certain 
orders. 

The compact finite difference (FD) method for solving Poisson equations has been well 
studied since 1984 [1,2,3,4,5,6,7,8,9,10,11,12,13]. Specifically, two-dimensional (2D) 
and 3D Poisson equations can be solved by high-order compact FD methods [1,2,3, 4,5, 
6]. These schemes are called “compact” since they only use minimum grid points to achieve 
fourth-order accuracy explicitly in the discretization formulas. Moreover, there has been a 
renewed interest in combining high-order compact scheme with multigrid method to solve 
Poisson equations. The classical multigrid method [14, 15, 16] combined with compact FD 
method for solving 2D and 3D Poisson equations has been conducted in [7,8,9, 10,11,12, 
13,17]. For example, Wang and Zhang [11] proposed a Richardson extrapolation for the 
numerical solutions from the two-level grids together with an operator based interpolation 
iterative strategy to achieve sixth-order accuracy by using the classical multigrid method 
and the fourth-order compact FD scheme. Ge [13] developed a fourth-order compact FD 
method with the classical multigrid method to solve the 3D Poisson equation using unequal 
meshsizes in different coordinate directions. Dehghan et al. [17] solved the ID, 2D and 3D 
Poisson equations with both second-order and fourth-order compact FD methods based on 
a new two-grid multigrid method. Besides Poisson equation, the classical multigrid method 
has been applied to many problems, including the biharmonic equation [18], the convection- 
diffusion equation [19,20,21] and so on. 

Cascadic multigrid (CMG) method proposed by Deuflhard and Bornemann in [22] is a 
variant of the multigrid without any coarse grid correction steps, where instead of starting 
from the finest grid, the solution is first computed on the coarsest grid and the recursively 
interpolated and relaxed on finer grids. Bornemann and Deuflhard [22] showed that it is an 
optimal iteration method with respect to the energy norm. Since the 1990s, the method has 
been frequently used to solve the elliptic equation with the finite element (FE) discretization 
because of its high efficiency and simplicity [24,7,26,?,?,?,?,?,?,?,?]. In 2007, Shi et al. [35] 
proposed an economical cascadic multigrid method using the different criteria for choosing 
the smoothing steps on each level of grid. Later, based on a new Richardson extrapolation 
formula for the linear FE solution, an extrapolation cascadic multigrid (EXCMG) method 
was first proposed by Chen et al. [36, 37] to solve 2D Poisson equation with the linear FE 
discretization. For the EXCMG method, in order to obtain a better initial guess of the iter¬ 
ative solution on the next finer grid, numerical solutions on the two-level of grids (current 
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and previous grids) are needed (whereas only one-level of numerical solution is needed in 
the CMG method). The EXCMG algorithm has been successfully applied to non-smooth 
problems [38], linear parabolic problems [39], and the simulation of the electric field with 
a point singularity arising in geophysical exploration [40,41], However, to our best knowl¬ 
edge, the EXCMG algorithm has mainly been used for solving the 2D elliptic problems with 
the linear FE discretization in existing literature. But it is of more importance to solve the 
3D problems efficiently and accurately arising in many engineering areas, such as geophys¬ 
ical exploration [42]. Since the construction process of the higher-order (at least fifth-order) 
approximation to the fourth-order compact FD solution on the next finer grid has to be dif¬ 
ferent from the construction process of the third-order approximation to the second-order FE 
solution, it will be nontrivial to extend the EXCMG method from 2D problems with second- 
order FE discretization to 3D problems with fourth-order compact FD discretization. 

In this paper, we will propose an EXCMG method combined with the fourth-order com¬ 
pact difference scheme to solve the Dirichlet boundary value problem of the 3D Poisson 
equation (1) in rectangular domains. In our approach, the computational domain is dis¬ 
cretized by regular grids, and a 19-point fourth-order compact difference scheme is used 
to discretize the 3D Poisson equation with unequal meshsizes in different directions. By 
combining the Richardson extrapolation and tri-quartic Lagrange interpolation for the nu¬ 
merical solutions from two-level of grids (current and previous grids), we are able to obtain 
a much better initial guess of the iterative solution on the next finer grid than one obtained by 
using linear interpolation in CMG method. Then, the resulting large linear system is solved 
by the conjugate gradient (CG) solver using the above obtained initial guess. Additionally, a 
tolerance related to relative residual is introduced in the CG solver in order to obtain conve¬ 
niently the numerical solution with the desired accuracy. Moreover, when the exact solution 
is sufficiently smooth, a simple method based on the midpoint extrapolation formula can 
be used to obtain cheaply and directly a sixth-order accurate solution on the entire fine grid 
from two fourth-order FD solutions on two different scale grids (current and previous grids). 
And a fourth-order compact FD scheme can be used to compute the gradient of the solution 
by solving a series of tridiagonal linear systems. Finally, our method has been used to solve 
3D Poisson equations with more than 16 million unknowns in about 10 seconds on a desk¬ 
top with 16GB RAM installed, which is much more efficient than the classical multigrid 
methods. 

The rest of the paper is organized as follows: section 2 gives the description of the com¬ 
pact FD discretization for the 3D Poisson equation. Section 3 reviews the classical V-cycle 
and W-cycle multigrid methods. In section 4, we first derive some sixth-order extrapolation 
formulas, and then develop a new EXCMG method to solve 3D Poisson equation. Section 5 
presents the numerical results to demonstrate the high efficiency and accuracy of the pro¬ 
posed method. And conclusions are given in the final section. 


2 Compact difference scheme 


We consider a cubic domain Q = [0, L r ] x [0, Ly\ x [0, L ; ], and discretize the do¬ 
main with unequal meshsizes h x , h y and h z in the x, y and z coordinate directions, respec¬ 
tively. Let N x = Lx/h x , N y = L y /h y , N z = L z /h z be the numbers of uniform intervals 
along the x, y and z directions. The grid points are ( X\,yj,Zk ), with x, = ih x ,yj = jh y and 
Zk = kh z , i = 0,1, ■ • • , N x , j = 0,1, • ■ ■ , N y and k = 0,1, - ■ ■ , N z . The quantity represents 
the numerical solution at (x,-, yj, Zk )■ 
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Then the value on the boundary points Uijjti = 0, N x or j = 0, N y or k = 0, N z ) can 
be evaluated directly from the Dirichlet boundary condition. For internal grid points (i = 
1, • • • , N x — 1, j = 1, • • • , N y - 1, k = 1, • • • , N z — 1), the 19-point fourth-order compact 
difference scheme with unequal-meshsize for 3D Poisson equation was derived in [6,13]: 
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~ 2 fi+l,j,k + fi-l,j,k + fi,j+l,k + fi,j-\,k + fi,j,k -1 + fi,j,k+ 1)* 


(*<i+lj+l,* + Ui+l,j-l,k + Ui-l,j+l,k + 


(2) 


Let li = ma x[h x ,h y ,h z ), throughout this paper, we denote u/, to be the FD solution of 
(2) with mesh sizes h x ,h y ,h z , while use u /,/2 to denote the FD solution of (2) when mesh 
sizes are h x / 2, h y /2, h-J 2. Then the difference scheme (2) can be expressed in the following 
matrix form: 

A h u h = //„ (3) 

where A/, is a sparse positive definite matrix, and //, denotes the right hand-side vector of (2) 
with mesh sizes li x , h y and h z . 


3 Classical multigrid methods 

The multigrid method is based on the idea that classical relaxation methods strongly damp 
the oscillatory error components, but converge slowly for smooth error components [15, 
16]. Flence, after a few relaxation sweeps, we compute the smooth residual of the current 
approximation V;, (with mesh sizes h x , h y , /;-) and transfer it to a coarser grid Q 2 h (with mesh 
sizes 2h x ,2h y ,2h z ) by a restriction operation, where the errors become more oscillatory. 
Solving the residual equation on the coarse grid Q 2 h, interpolating the correction back to the 
fine grid Qi,, and adding it to the fine-grid current approximation v/, yields to the two-grid 
correction method. Since the coarse-grid problem is not much different from the original 
problem, we can perform a few, say y, two-grid iteration steps (see Fig. 1) to the residual 
equation on the coarse grid, which means relaxing there and then moving to Qy, (with 
mesh sizes 4 h x , 4h y , 4h z ) for the correction step. We can repeat this process on successively 
coarser grids until a direct solution of the residual equation is possible. Then the corrections 
are interpolated back to finer grids until the process reaches the finest grid Qj, (with mesh 
sizes h x , h y , h z ) and the fine-grid approximate solution is corrected. 

Usually, the cases y = 1 and y = 2 are particularly interesting. We refer to the case 
y = 1 as V-cycle and to y = 2 as W-cycle. The number y is also called cycle index. A V- 
cycle multigrid method is obtained when the V-cycle is repeated until a stopping criterion is 
satisfied on the finest grid. We refer to a V-cycle (W-cycle) with v\ relaxation sweeps before 
the correction step and v 2 relaxation sweeps after the correction step as a V(vi, v 2 )-cycle 
(W(vi, v 2 )-cycle). 
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Fig. 1 The four-level structure of the V-cycle, W-cycle, FMG, CMG and EXCMG methods. In the diagram, 
• denotes pre-smoothing, o denotes post-smoothing, ft denotes prolongation (usually defined by linear inter¬ 
polation), i denotes restriction, ft denotes extrapolation and high-order interpolation, and ■ denotes direct 
solver. 


4 Extrapolation cascadic multigrid methods 


The CMG method proposed by Deuflhard and Bornemann in [22] is a variant of full multi¬ 
grid (FMG) method without any coarse grid correction steps but with an a posteriori con¬ 
trol of the number of smoothing iterations (see Fig. 1). It has been shown that the CMG 
method has optimal computational complexity for both conforming and nonconforming el¬ 
ements with CG as a smoother. Since the 1990s, the CMG method has received quite a 
bit of attention from researchers because of its high efficiency and simplicity [23,24,7,26, 
?,?,?,?,?,?,?,?]. 

In 2008, by using Richardson extrapolation and bilinear quadratic interpolation for the 
FE solutions on two-level of grids (current and previous grids) to obtain an extremely ac¬ 
curate initial guess of the iterative solution on the next finer grid, Chen et al. [36] proposed 
an extrapolation cascadic multigrid (EXCMG) method to solve 2D elliptic boundary value 
problems. It has been shown in [37] that the EXCMG method is much more efficient than the 
CMG method, which simply uses the linear interpolation for the FE solution on the current 
grid to provide an initial guess of the iterative solution on the next finer grid. Recently, we 
improved and generalized the EXCMG method to solve large linear systems resulting from 
FE discretization of 3D elliptic problems, compared it with the classical multigrid methods, 
and further presented the reason why EXCMG algorithms are highly efficient [44]. How¬ 
ever, to our best knowledge, CMG and EXCMG are mainly used for linear FE method in 
existing literature, and it will be interesting to extend the EXCMG method to the field of 
high-order FD method. 
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4.1 EXCMG algorithm combined with compact difference scheme 

The key ingredients of the EXCMG method are extrapolation and high-order interpolation 
(see Fig. 1), which can produce a much better initial guess of the iterative solution on the 
next finer grid than one obtained by using linear interpolation in CMG method. 

In this subsection, we will propose a new EXCMG method combined with fourth-order 
compact difference scheme for solving the Dirichlet boundary value problem of the 3D 
Poisson equation, which is stated in the following algorithm. 


Algorithm 1 New EXCMG method: u h ) <= EXCMG(A/,, //,, L, e) 

1: uu <= DSOLVE(A„m /7 = [a ) > uh is FD solution of (3) with mesh sizes //,, H ,, H-. 

2: uh /2 <= It SOI .V \[( A njiUj i/ 2 = fun) > »a /2 is FD solution of (3) with mesh sizes H x j 2, H y /1, HJ2. 

3: h x = H x /2, h y = H y /2, h z = H-J2 

4: for i = 1 to L do 

5: h x = h x /2, h y = h y /2, h z = hJ2 

6: Wk = EXPfi„it e (u 2 h, t> tv/, is a fifth-order approximation of the actual numerical solu¬ 

tion Uj,. and it serves as the initial guess for the CG solver on the next finer grid. 

7 : while \\A h u h - /;,|| 2 > e ■ ||/,,|| 2 do 

8 : ui, <= CG(A h ,u h ,f h ) 

9: end while 

10: = HXP,, U 2 h) > Optional step. «/, is a sixth-order approximation solution for 

sufficiently smooth 11 . 

11: end for 


In Algorithm 1, the coarsest grid has the mesh sizes H x ,H y ,H z , the positive integer L 
is the total number of grids except first two embedded grids, which indicates that the mesh 
sizes of the finest grid are DSOLVE is a direct solver used on the first two 

coarse grids (see line 1-2 in Algorithm 1). Procedure ¥XPf inite {u 2 h, « 4 /i) denotes a fifth-order 
approximation to the actual compact FD solution m/, obtained by Richardson extrapolation 
and tri-quartic Lagrange interpolation from the numerical solutions wj/, and u 4 /,. And there is 
an optional step in the above algorithm (see line 10 in Algorithm 1), where EXP frw (M/,, U 2 h) 
denotes a higher-order approximation solution on entire fine grid with mesh size h from 
two fourth-order FD solutions m/, and i< 2 /,. This optional step is used to increase the order of 
solution accuracy from fourth order to sixth order (see Table 1-10 in section 5 for details) 
when the exact solution u of elliptic equation ( 1 ) is sufficiently smooth. 

The detailed procedures of extrapolation and tri-quartic Lagrange interpolation are de¬ 
scribed in the next two subsections 4.2 and 4.3. The differences between our new EXCMG 
method and existing EXCMG method [36,37] are listed as follows: 

(1) In our new EXCMG method, a fourth-order compact difference scheme, rather than the 
second-order linear FE method, is employed to discretize the 3D Poisson equation. 

(2) Instead of performing a fixed number of smoothing iterations as used in the existing 
EXCMG method [36,37], a relative residual tolerance e is introduced for the smoother 
in our EXCMG method (see line 7 in Algorithm 1), which enables us to conveniently 
obtain the numerical solution with the desired accuracy. 

(3) In the existing EXCMG literature [36,37], a third-order approximation to the second- 
order FE solution is constructed to serve as the initial guess for the iterative solver on the 
next finer grid, and the construction of the third-order approximation to the second-order 
FE solution is done at every single coarse hexahedral element. However, in our new 
EXCMG method, a fifth-order approximation to the fourth-order FD solution, obtained 
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through the Richardson extrapolation and tri-quartic Lagrange interpolation, is used as 
the initial guess for the iterative solver. In addition, the tri-quartic interpolation should 
be done for every cell which contains eight neighboring coarse hexahedral elements as 
shown in Fig. 3, rather than every single coarse hexahedral element. 


4.2 Extrapolation and quartic interpolation: ID case 

The extrapolation method is an efficient procedure for increasing the solution accuracy of 
many problems in numerical analysis. Marchuk and Shaidurov [45] systematically studied 
its application in the FD method in 1983. Since then, this technique has been well demon¬ 
strated in the framework of the FD and FE methods [46,47,48,21,49,50,51,52,53,54], 

In this and next subsections, we assume that the exact solution u is sufficiently smooth, 
and we will formally explain how to use extrapolation and quartic interpolation techniques 
to obtain the fifth-order approximation Wf, of the fourth-order FD solution on the next finer 
grid, which can be regarded as another important application of the extrapolation method. 
In addition, we will also show how to construct the enhanced sixth-order accurate numerical 
solution Uh for the problem (1). 

4.2.1 Extrapolation for the true solution 

For simplicity, we first consider the three-levels of embedded grids Zfi = 0,1,2) with mesh 
sizes lii = ho/2' in one dimension. Suppose u e H 6 (Q), from theorem 4.1 in [43] (taking 
m = 2, s = 6) and by using the result that H 2 (Q) can be continuously embedded into 
we can get that the error ||e'||oo should be 0(h 4 ), where e‘ = u' — u is the error of the fourth- 
order compact FD solution u' with mesh size hi. Now we further assume that the truncation 
error at node x k has the form 


e\x k ) = A(x k )hf + O(hf), (4) 

where A(x) is a suitably smooth function independent of h, . The truncation error expansion 
(4) will be verified by numerical results in section 5. 

It is well known that the traditional extrapolation is possible only at coarse grid points, 
where at least two approximations, corresponding to different mesh size, are known. From 
eq. (4), we easily obtain the Richardson extrapolation formula at coarse grid points 

16 u! — m? . 

u\ :=-—-= u(x k ) + O(h 0 ), k = j,j + 1, (5) 

which is a sixth-order approximation to the true solution at the coarse grid points. 

In fact, by using the linear interpolation formula, one can also obtain a sixth-order ac¬ 
curate approximation at the fine grid point xj+ 1 / 2 - Setting i = 0 and i = 1 in eq. (4) and then 
subtracting each other, we have 

Mxic) = tttt ( M ° - u k) + OQil), k = j, j + 1 . ( 6 ) 

15 h 0 

From the error estimate of the linear interpolation 

1 

Mx j+ 1 / 2 ) = -(A(xj) + A(x j+ 1 )) + O(h~ 0 ), 


(7) 
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Fig. 2 Three embedded grids for two neighboring coarse elements in ID. 


and substituting eq. ( 6 ) into eq. (7) , we get 

Mxj+1/2) = u°j - “]) + ^(w° + 1 - “]+t) + 0(/!5). (8) 

Since 

“]+ 1/2 = “(* 7 + 1 / 2 ) + -j^U* 7 + 1/2 )/lo + (9) 

by using ( 8 ), we obtain the following midpoint extrapolation formula: 

“}+i/2 : = «] + i/2 + 3 q(“) - “/ + “l+i - “/+1) = “(*7+1/2) + 0(/1 q), (10 ) 

which is a sixth-order approximation to the true solution at the fine grid point X/+1/2- 


4.2.2 Extrapolation for the FD solution 

In this part, we will explain, given the fourth-order FD solutions u° and u 1 , how to use 
the extrapolation and high-order interpolation techniques to construct a fifth-order (to be 
illustrated in subsection 4.4) approximation w 2 to the FD solution n 2 . 

Adding one midpoint and two four equal division points, the coarse mesh element 
(xj, Xj+ 1 ) is uniformly refined into four elements of fine mesh Z 2 as shown in Fig. 2. Assume 
there exists a constant c such that 

CM 1 + (1 - C)u° = U 2 + 0(/1q). (11) 

Here, we aim to use a linear combination of u° and if to approximate the FD solution if up 
to sixth-order accuracy. Substituting the asymptotic error expansion (4) into (11), we obtain 
c = 17/16 and an extrapolation formula 

17 if — iP 

w 2 :=-= u 2 + O(h 6 0 ), k = j,j+ 1, (12) 

at nodes xj and Xj+ \. To derive the extrapolation formula at midpoint Xj+ 1 / 2 , eq. (4) leads to 

“J+1/2 = “j+1/2 ~ 25ftA( x j+il2)hQ + OQPf). 


( 13 ) 
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Zo 


Zi 



Fig. 3 Three embedded grids on one interpolation cell which contains eight neighboring coarse hexahedral 
elements. 


Substituting eq. (8) into eq. (13), we have the following sixth-order extrapolation formula at 
the midpoint Xj+ 1 / 2 , 

w f+i /2 := m /+i /2 + 32 *-“/ - u i + u i + 1 - u j+ = u j+ 1/2 + O(h 6 0 ). (14) 

Sixth-order extrapolation formulas (12) and (14) can be efficiently applied to each coarse- 
grid element (xj, Xj+ 1 ). 

Once the five approximated values vv 2 ,vv 2 +1( , 2 ,vv 2 +1 ,vv 2 +3( , 2 an d w /+2 arc obtained on the 
two neighboring coarse elements, we can get the following four equal division point extrap¬ 
olation formulas by using the quartic interpolation 


W /+l/4 

:= 1^8 '( 35Vl ; + 140VV 1+1/2 “ 7 ° W i+l + 28w f+3/2 - Sw ) + 2 )’ 

(15) 

W )+ 3/4 

:= 128 ~ ^ + 60w i +1 /2 + 90w }+i _ 20wl j+ 3/2 + 3m ; J+2)> 

(16) 

W )+ 5/4 

:= 1^8 '( _ 5wZ J + 28w 7 +i /2 - 1Qw ) + \ + U0w2 j+ 3/2 + 35w2 j+ 2 ), 

(17) 

W J+7/4 

: = ( 3w i _ 2Qw ) + 1/2 + 90w2 j + 1 + 60w f + 3/2 - 5w2 j + i)- 

(18) 


Until now, we have obtained a high-order approximation w 2 to the FD solution u 2 , which 
can be used as the initial guess of the iterative solution on the fine mesh Z 2 . 



4.3 Extrapolation and quartic interpolation: 3D case 

In this subsection, we explain how to obtain a fifth-order accurate approximation w 2 to 
the fourth-order FD solution u 2 , and a sixth-order accurate approximate solution u l to the 
problem (1) for embedded hexahedral grids as shown in Fig. 3. 

Taking every interpolation cell which consists of eight neighboring coarse hexahedral 
elements (see Fig. 3) into account, the construction processes of the approximation w 2 are 
as follows: 

Corner Nodes (such as 1, 3, 51, 53): The approximate values at 27 corner nodes *•’ on such 
interpolation cell can be obtained by using the extrapolation formula (12). 
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Midpoints of edges (such as 2, 6, 26, 28): The approximate values at these 54 midpoints 
on such interpolation cell can be obtained by using the midpoint extrapolation for¬ 
mula (14) in x-direction, y-dircction or z-direction. 

Centers of faces (such as 27, 31, 107, 109): Since the center of each face on such interpo¬ 
lation cell can be viewed as the midpoint of two face diagonals, using the midpoint ex¬ 
trapolation formula (14) we can obtain two approximate values, and take the arithmetic 
mean as the approximation at these 36 midpoints . 

Centers of coarse hexahedral elements (such as 32, 42, 82, 92): Since the center of each coarse 
hexahedral element on such interpolation cell can be viewed as the midpoint of four 
space diagonals, again using the midpoint extrapolation formula (14) we can obtain four 
approximate values, and take the arithmetic mean as the approximation at these 8 mid¬ 
points . 

Other fine grid points: The approximate values of remaining 604(9 3 - 5 3 ) fine grid points 
on such the interpolation cell can be obtained by using tri-quartic Lagrange interpolation 
with the known 125-node (27 corner nodes, 54 midpoints of edges, 36 centers of faces 
and 8 centers of coarse hexahedral elements) values. 


is 


The tri-quartic Lagrange interpolation function in terms of natural coordinates (£, r/, £) 

125 

w 2 {f,r 1 ,0 = y j N m (f,i h Owl, (19) 

m= 1 


where the shape functions N m can be written as follows 


Nm(.Z,ri,0 = /?(f)/-(??)4(a 


( 20 ) 


where l*(x) (0 < i < 4) is the Lagrange fundamental polynomials of degree 4, defined as 


?©= n 

k=0,k*i 


£i - ft ’ 


( 21 ) 


and (£,-, rjj, £ k ) is the natural coordinate of node m (1 <m< 125). 

When constructing the sixth-order accurate solution u l based on two fourth-order accu¬ 
rate solutions u° and u 1 , the Richardson extrapolation formula (5) can be directly used for 
coarse grid points, while the sixth-order midpoint extrapolation formula (10) can be directly 
used for all other fine grid points, which is similar to the process (excluding the tri-quartic 
interpolation) of constructing the approximation ir 2 described as above. 


Remark 1 Since the compact FD solution »/, of (2) is a fourth-order approximation of the 
exact solution u, in order to get a quite good initial guess w/, for the CG solver, a tri-quartic 
Lagrange interpolation method is employed in this paper so that a fifth-order approximation 
of w/, to ui, is achieved. Moreover, the relative effect of how w/, approximates m/, becomes 
better when mesh is refined, thus, the number of iterations will be reduced most signifi¬ 
cantly on the finest grid, which is particularly important for solving large linear systems 
and can greatly reduce the computational cost. We note that the tri-quadratic interpolation 
used in [44] produces a third-order approximation to the second-order FE solution, and the 
tri-quadratic interpolation is accurate enough in that case. Flowever, when u/, is obtained 
from the fourth-order compact FD method as shown in this paper, the tri-quadratic interpo¬ 
lation can not provide a sufficiently accurate initial guess w/,, the relative effect of how w/, 
approximates «/, will become worse when mesh is refined. 
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Remark 2 Tri-quartic Lagrange interpolation defined by eq. (19) is a local operation defin¬ 
ing on each interpolation cell containing eight neighbouring coarse elements. In fact, eq. 
(19) defines a same (604 x 125) interpolation matrix on every interpolation cell, thus the ap¬ 
proximate values of remaining 604(9 3 - 5 3 ) fine-grid points on every interpolation cell can 
be obtained by multiplying the (604 x 125) interpolation matrix with the vector consisting of 
125 known values on such interpolation cell. Therefore, the fifth-order approximation of FD 
solution w/, on the entire domain can be obtained very effectively by applying the extrapo¬ 
lation formulas (12) and (14) to the 125 nodes mentioned above, and running the tri-quartic 
Lagrange interpolation (19) based on such 125 known values for every interpolation cell in 
the entire domain. 


4.4 The error analysis of initial guess w 2 

Let e = w 2 - u 1 be the difference between the initial guess w 2 and the FD solution u 2 . 
Assume that e has continuous derivatives up to order 5 on interval [xj,Xj + 2 ]. From (12) and 
(14) we obtain the equation 

e(x k ) = O(h 6 0 ), k = j, j + 1/2, j + 1, j + 3/2, j + 2. (22) 

From polynomial interpolation theory, the error of quartic interpolation I\ f can be repre¬ 
sented as 


Ra(x) = e-he= —e (5) (£)(x - xj)(x - x j+ U2 )(x - x j+ i)(x - x j+3/2 )(x - x j+2 ), 
where £ e (xj, Xj +2 ) depends on x. Especially at four equal division points we have 


Ra(xj+ 1/4) - 


Ra(Xj+3/a) - - 


21 i_ 

8 x 4 5 


Ihl 


6 tfl) = si92 e ? Kx j+ x) + °iK\ 

S x45' ra * ,= -ti> « <5> <^ >♦»<'*• 


3 hi 


and 


Ra(Xj+5/a) - 


3 hi 


8 x 4 5 
Ihl 


3 h 5 

0 e (5) (f 3 ) = Fi71T e<5) ( x ;+i) + °( /z o) ~ ^( xj ^/ a ). 


8192 
7 h 5 n 


Ra(x j+ 7 / 4 ) = -^ 5 e(5> (&) = -gj^e ( 5 ) Uf + t) + o(h 5 0 ) « -R A (x j+ 1 / 4 ). 
It follows from eqs. (22) and (24)-(27) that 

e{x k ) = I A e(x k ) + R A (x k ) = O(h 5 0 ), k = j + 1/4, j + 3/4, j + 5/4, j + 7/4, 


(23) 

(24) 

(25) 

(26) 

(27) 

(28) 


which means that the initial guess w 2 obtained by extrapolation and quartic interpolation is 
a fifth-order accurate approximation to the FD solution u 2 . 

The above error analysis can be directly extended to 3D case (see numerical verification 
in Section 5: the last columns in Table 1-11). In addition, eqs. (26) and (27) imply that the 
initial error e{x) forms a high-frequency oscillation in the entire domain, however, it can be 
smoothed out after a few CG iterations (see Fig. 4 for details). 
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5 Numerical experiments 


In this section, in order to illustrate the efficiency of the new EXCMG method comparing to 
the classical V-cycle and W-cycle multigrid methods with the Gauss-Seidel relaxation and 
the CG relaxation, we present the numerical results for six examples with smooth and finite 
regular solutions using the proposed method. Our code is written in Fortran 90 with double 
precision arithmetic, and compiled with Intel Visual Fortran Compiler XE 12.1 under 64- 
bit Windows 7. All programs are carried out on a personal desktop equipped with Intel(R) 
Core(TM) i7-4790K CPU (4.00 GHz) and 16GB RAM. 

The order of convergence of the method is computed by 


order = log 2 


IN - u\\ 

iN/2- wir 


(29) 


where || • || denotes some norm (for instance, Zr-norm or L“-norm) and u is the true solution. 


5.1 Numerical accuracy 

Example 1 The test Problem 1 can be written as 

d 2 u d 2 u d 2 u to s 

^2 + fyi + = e sin(xy)(l - jc 2 -y 2 ), in Q = [0,1] , 

where the boundary conditions are 

u{ 0 , y, z) = u(x, 0, z) = 0, u{ 1 , y, z) = e z sin(y), u(x, 1 , z) = e z sin(x), 

and 


(30) 


u{x,y, 0 ) = sin(xy), u(x,y, 1 ) = e sin(xy). 
The analytic solution of eq. (30) is 

u(x,y,z ) = e z sin(jcy). 


which is a sufficiently smooth function. 


Using 7 embedded grids with the coarsest grid 4 x 4 x 4, we present the numerical 
results for Problem 1 obtained by the new EXCMG method with e = 10 -14 in Table 1-2. 
Table 1 lists the Z. 2 -error of the compact FD solution the Lr -error of the gradient of the 
FD solution Vm;„ the L 2 -error of the extrapolated solution S/,, the i 2 -norm of the difference 
between the initial guess w/, and the FD solution «/,, and corresponding convergence rates. 
Table 2 gives all errors and convergence rates in L“-norm. Since a direct solver is used for 
the first two coarse levels of grids, we only list the results starting from the third level of 
grid 16 x 16 x 16. 

Here we explain how to numerically compute the gradient V it/, after we obtain the FD 
solution it/,. First, we use the following fourth-order, one-sided, FD approximation of the 
partial derivative u x on the boundary grid points, 


25 4 3 4 1 , , 

(^jc)o ,j,k — \2h~ U °’j ,k ~h ~ — ~fi ~^3 ,j,k — 4 ^ ^4 ,j,k> j — 0, • • • , Ny, k — 0, • • • , N z , 


25 


3 


1 


{Ux)N x ,j,k - 12/i U Nx,j,k - J^ u N x -\,j,k + -^-UN x -2,j,k ~ ,jjv for j - 0, ■ ■ ■ , N y ,k - 0. 


,N Z . 
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Table 1 Errors and convergence rates with e = 10 14 in Zr-norm for Example 1. 


mesh 

I|h/, - «||2 

order 

l|V(«* - u )|| 2 

order 

Pa - w||2 

order 

llw* - «aI|2 

order 

16 x 16 x 16 

1.67(—08) 


1.08(—06) 


1.381-09) 


4.361-07) 


32x 32x 32 

1.09(—09) 

3.93 

4.741-08) 

4.51 

2.401-11) 

5.84 

1.291-08) 

5.08 

64 x 64 x 64 

7.00(—11) 

3.97 

2.121-09) 

4.48 

3.911-13) 

5.94 

3.941-10) 

5.04 

128 x 128 x 128 

4.42(—12) 

3.98 

9.701-11) 

4.45 

6.221-15) 

5.98 

1.221-11) 

5.02 

256 x 256 x 256 

2.82(—13) 

3.97 

4.621-12) 

4.39 

5.681-15) 

0.13 

3.791-13) 

5.01 

Table 2 Errors and 

convergence rates with e 

= 10- 14 in L°°- 

■norm for Example 1. 




mesh 

I|h/i - »IU 

order 

l|V(«A - u )Hoc 

, order 

Pa - «l|oo 

order 

\\Wh - MAlloo 

order 

16 x 16 x 16 

5.47(—08) 


8.251-06) 


9.12(-09) 


3.851-06) 


32 x 32 x 32 

3.43(—09) 

4.00 

5.181-07) 

3.99 

1.771-10) 

5.68 

1.311-07) 

4.88 

64 x 64 x 64 

2.15(—10) 

4.00 

3.241-08) 

4.00 

3.201-12) 

5.79 

4.191-09) 

4.96 

128 x 128 x 128 

1.34(—11) 

4.00 

2.031-09) 

4.00 

5.421-14) 

5.88 

1.321-10) 

4.99 

256 x 256 x 256 

8.49(—13) 

4.00 

1.271-10) 

4.00 

4.001-14) 

0.44 

4.151-12) 

4.99 


Then we can obtain (u x \jj„(i = 1, • ■ ■ , N x — 1) on the internal grid points by solving the 
following linear system resulting from the fourth-order compact FD scheme [55], 


1 

6 


4 K 

1 J,k 3 " ~^(Mx)i,j,k 3 " 


2 hx 


for j = (),■•• ,N y ,k = (),■■■ ,N Z . 


The above ID tridiagonal system can be solved fast by the Thomas algorithm. Clearly, we 
can get u y and u z from similar procedures. Then, V«/, can be obtained efficiently. 

As we can see from table 1-2 that initial guess Wh is a fifth-order approximation to the 
FD solution Uh, which validates our theoretical analysis in section 4.4, and the FD solution 
u/j achieves the full fourth-order accuracy. The numerical gradient Vnj, is also a fourth-order 
approximation to the exact gradient V« in both the Zr-norm and L“-norm, while the extrap¬ 
olated solution 5/, converges with sixth-order accuracy on all grids except the finest grid. 
This is due to the fact that the extrapolated solution zJ/, is obtained from two fourth-order 
FD solutions «/, and u 2 /,, these two solutions must be extremely accurate in order to obtain 
a sixth-order accurate solution u/,. As the grid becomes finer, the relative residual tolerance 
needs to be smaller. Thus, the extrapolated solution »/, starts to lose convergence order when 
the grid is fine enough since a uniform tolerance is used in our EXCMG algorithm. And in 
this example, on the finest mesh 256 x 256 x 256, the maximum error between the extrap¬ 
olated solution S/, and the exact solution u already reaches O(10 -14 ), which is almost the 
machine accuracy, although the method does not achieve the full sixth-order on the finest 
grid. Additionally, we can see that the numerical results confirm with the asymptotic error 
expansion (4). 


Example 2 The test Problem 2 can be written as 


d 2 u d 2 u d 2 u 
dx 2 dy 2 dz 2 


in Q = [0,1]\ 


( 31 ) 


where the boundary conditions are 


u( 0, y, z) = e y sin( V2 z), n( x, 0, z) = e x sin( V2 z), u(x, y, 0) = 0, 


and 


x+> sin( V2z), u(x,y, 1) = e x+y sin( V2). 


u{l,y,z) = e 1+y sin(V2z), u(x,l,z) = e 
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Table 3 Errors and convergence rates with 6=10 14 in L 2 -norm for Example 2. 


mesh 

I|h/i - u\h 

order 

l|V(«A - U) || 2 

order 

II- «ll2 

order 

IK - u h lb 

order 

16 x 16 x 16 

4.26(—08) 


150(—05) 


2.28(-09) 


6.21(-06) 


32 x 32 x 32 

2.79(—09) 

3.94 

663(—07) 

4.50 

3.886-11) 

5.88 

1.956-07) 

5.00 

64 x 64 x 64 

1.78(—10) 

3.97 

295(—08) 

4.49 

6.28(-13) 

5.95 

6.106-09) 

5.00 

128 x 128 x 128 

1.13(—11) 

3.98 

132(—09) 

4.48 

1.026-14) 

5.94 

1.916-10) 

5.00 

256 x 256 x 256 

7.32(—13) 

3.94 

596(—11) 

4.47 

3.116-14) 

-1.60 

5.976-12) 

5.00 

Table 4 Errors and 

convergence rates with e 

= 10^ 14 in L°°- 

■norm for Example 2. 




mesh 

I|h/i - «IU 

order 

l|V(«A - «)||oc 

, order 

II u h ~ W||oo 

order 

IK - «;,||oo 

order 

16 x 16 x 16 

1.16C-07) 


1.26(-4) 


1.016-08) 


3.396-05) 


32 x 32 x 32 

7.226—09) 

4.00 

7.95(—6) 

3.99 

1.866-10) 

5.76 

1.116-06) 

4.94 

64 x 64 x 64 

4.52(—10) 

4.00 

4.98(—7) 

4.00 

3.176-12) 

5.87 

3.546-08) 

4.96 

128 x 128 x 128 

2.83(— 11) 

4.00 

3.116—8) 

4.00 

5.956-14) 

5.74 

1.126-09) 

4.99 

256 x 256 x 256 

1.81(—12) 

3.96 

1.951-9) 

4.00 

1.076-13) 

-0.85 

3.516-11) 

4.99 


The analytic solution of eq. (3 1 ) is 

u = e x+y sin( V2z), 

which is a harmonic function and has arbitrary order smooth derivatives. 

Again, we use 7 embedded grids with the coarsest grid 4x4x4, and the corresponding 
numerical results obtained by the EXCMG method with e = 10“ 14 are listed in table 3 
and 4. Once again, initial guess w /, is a fifth-order approximation of the FD solution it/,, the 
FD solution is fourth-order accurate, and the numerical gradient V«/, is also a fourth-order 
approximation to the exact gradient V//, while the extrapolated solution it/, converges to exact 
solution u with sixth-order but starts to lose accuracy on the finest grid 256 x 256 x 256. 
Additionally, the maximum error between the extrapolated solution it/, and the exact solution 
u is less than 6.0 x 10 -14 , which means that the solution iti, is already accurate enough, and 
we don’t need to further reduce the error tolerance. 


Example 3 The test Problem 3 can be written as 


d 2 u d 2 u d 2 u o 

37 + 57 + 3? = /(w> - mG = [<UI - 

u = g(x, y, z), on dQ, 


where / and g are determined from the exact solution 

u _ e -3((.v-0.5) 2 +(y-0.5) 2 +(z-0.5) 2 ) 


(32) 


which is a 3D Gaussian function. It varies rapidly near the point (0.5,0.5,0.5). 

Table 5 and 6 list the numerical results obtained by the EXCMG method with e = 10 _11 . 
One more time, one can see that initial guess wt, is a fifth-order approximation of the FD 
solution Uh, the FD solution ;</, is fourth-order accurate (although the convergent order is 
slightly reduced on the finest grid), and the numerical gradient Vu h is also a fourth-order 
approximation to the exact gradient Vh, while the extrapolated solution u/, is sixth-order 
accurate. Therefore, our EXCMG method is still very effective for the problem with very 
rapid variations. 
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Table 5 Errors and convergence rates with e = 10 11 in Zr-norm for Example 3. 


mesh 

II u h - «||2 

order 

l|V(«A - «)ll2 

order 

Pa - “Ik 

order 

\\Wh ~ Uhlk 

order 

16 x 16 x 16 

1.22(—06) 


4.08(-04) 


6.68(-08) 


2.29(-04) 


32 x 32 x 32 

7.90(—08) 

3.95 

1.29(—05) 

4.99 

1.15(—09) 

5.86 

5.80(-06) 

5.30 

64 x 64 x 64 

5.04(-09) 

3.97 

4.42(—07) 

4.86 

1.87(—11) 

5.94 

1.86(-07) 

4.96 

128 x 128 x 128 

3.19(—10) 

3.98 

1.73(—08) 

4.67 

3.05(—13) 

5.94 

5.86(-09) 

4.99 

256 x 256 x 256 

2.40(—11) 

3.73 

7.50(-10) 

4.53 

4.52(—12) 

-3.89 

1.84(—10) 

4.99 

Table 6 Errors and convergence rates with e 

= lO^ 11 in L°°- 

norm for Example 3. 




mesh 

pA ~ U\\„ 

order 

l|V (lih - u )IU 

order 

Pa - «||co 

order 

\\wi, - «aIU 

order 

16 x 16 x 16 

4.80(—06) 


1,03(—03) 


2.14(—07) 


1.13(—03) 


32 x 32 x 32 

2.97(—07) 

4.01 

4.35(-05) 

4.56 

4.07(—09) 

5.71 

2.50(-05) 

5.50 

64 x 64 x 64 

1.85(—08) 

4.00 

2.0K-06) 

4.43 

6.63(—11) 

5.94 

9.98(—07) 

4.65 

128 x 128 x 128 

1.16(—09) 

4.00 

1,03(—07) 

4.28 

1.04(—12) 

5.99 

3.02(-08) 

5.05 

256 x 256 x 256 

8.95(—11) 

3.69 

5.55(—09) 

4.22 

1.851-11) 

-4.15 

9.27(-10) 

5.03 

Table 7 Errors and convergence rates with e 

= 10 9 in L 2 -norm for Example 4. 




mesh 

II U h - «||2 

order 

I|V(«A - «)|| 2 

order 

Pa - “Ik 

order 

\\Wh - u h || 2 

order 

32 x 16 x 8 

3.58(—06) 


2.85(-4) 


4.55(-07) 


7.86(-05) 


64 x 32 x 16 

2.35(—07) 

3.93 

1.37(—5) 

4.38 

8.36(—09) 

5.76 

2.58(-06) 

4.93 

128 x 64 x 32 

1.51 (-08) 

3.96 

6.34(—7) 

4.43 

1.39(—10) 

5.92 

8.19(-08) 

4.98 

256 x 128 x 64 

9.51 (-10) 

3.98 

2.96(—8) 

4.42 

3.24(—12) 

5.42 

2.58(-09) 

4.99 

512 x 256 x 128 

5.73(—11) 

4.05 

1.67(—9) 

4.15 

9.12(—12) 

-1.49 

7.99(—11) 

5.01 


Example 4 The test Problem 4 can be written as 

d 2 u d 2 u d 2 u j n 

—r + —- + —- = -5.25 tt sin(27rjc) sin(ay) sin(-z), in Q = [0, l] 3 , (33) 

ax 1 dy z oz~ 2 

where the boundary conditions are 

u(0,y,z) = u(l,y,z) = u(x,0,z) = u(x, l,z) = u(x,y, 0) = 0 and u(x,y, 1) = sinilnx) sin(^ry). 
The analytic solution of eq. (33) is 


7T 

u(x,y,z) = sin(2fl , jt) sin(7ty) stn(— z), 

which changes fastest in the x direction, faster in the y direction and slowest in the z direc¬ 
tion. 


Since the solution has the fastest change in the x-direction and the slowest change in 
the z-direction, we use the coarsest grid 8 x 4 x 2 in the EXCMG algorithm. Table 7 and 8 
list the numerical data obtained by EXCMG method using a tolerance e = 10 -9 . Again, the 
initial guess w/, is a fifth-order approximation of the FD solution t</,, the FD solution m/, is 
fourth-order accurate, and the numerical gradient V«/, is also a fourth-order approximation 
to the exact gradient V«, while the extrapolated solution S/, achieves sixth-order accuracy 
but starts to lose accuracy on the finest grid since a uniform tolerance e = 10 -9 is used on 
each level of grid. The maximum error between the extrapolated solution it/, and the exact 
solution u already reaches O(10 _1 *) on the finest grid which is again quite accurate. 

Previous examples are results for the 3D Poisson equation where the exact solution is 
infinitely many times continuously differentiable. In the following examples, we will show 
the results using the new EXCMG method for the cases where the exact solutions have finite 
regularities. 
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Table 8 Errors and convergence rates with 6=10 9 in L°°-norm for Example 4. 


mesh 

\\Uh ~ »l|oo 

order 

l|V(«/, - «)||oo 

order 

With - w|U 

order 

II w h - U h ||oc 

order 

32 x 16 x 8 

1.10(—05) 


1.78(-3) 


2.741-06) 


3.581-04) 


64 x 32 x 16 

6.97(—07) 

3.98 

1.161-4) 

3.93 

6.231-08) 

5.46 

9.971-06) 

5.17 

128 x 64 x 32 

4.35(—08) 

4.00 

7.361-6) 

3.98 

1.181-09) 

5.72 

3.201-07) 

4.96 

256 x 128 x 64 

2.71(—09) 

4.01 

4.611-7) 

4.00 

1.921-11) 

5.94 

9.591-09) 

5.06 

512 x 256 x 128 

1,82(— 10) 

3.90 

3.741-8) 

3.63 

4.581-11) 

-1.26 

2.821-10) 

5.09 

Table 9 Errors and convergence rates with e = 

= 10 12 in L 2 -norm for Example 5. 




mesh 

II U h - «||2 

order 

I|V(«/, - n )||2 

order 

II Uh ~ »ll2 

order 

II W h - Uh\\l 

order 

16 x 16 x 16 

5,44(—08) 


1.22(-05) 


3.351-09) 


2.32(-06) 


32 x 32 x 32 

3.56(-09) 

3.94 

6.84(-07) 

4.15 

5.371-11) 

5.96 

1.291-07) 

4.18 

64 x 64 x 64 

2.27(—10) 

3.97 

3.251-08) 

4.39 

8.57(-13) 

5.97 

4.191-09) 

4.94 

128 x 128 x 128 

1.44(-11) 

3.98 

1.471-09) 

4.46 

1.341-14) 

6.00 

1.321-10) 

4.99 

256 x 256 x 256 

9.57(-13) 

3.91 

6.591-11) 

4.48 

1.181-13) 

-3.14 

4.121-12) 

5.00 

Table 10 Errors and 

convergence 

rates with e 

= 10~ 12 in L” 

-norm for Example 5. 




mesh 

ll«A - i'llco 

order 

l|V(«;, - «)||oo 

order 

II “h ~ "lice 

order 

II Wh - UhWoa 

order 

16 x 16 x 16 

1.131-07) 


8.261-5) 


1.181-08) 


2.531-05) 


32 x 32 x 32 

7.161-09) 

3.98 

5.831-6) 

3.82 

2.081-10) 

5.83 

6.961-07) 

5.18 

64 x 64 x 64 

4.481-10) 

4.00 

3.761-7) 

3.96 

3.32(-12) 

5.97 

2.481-08) 

4.81 

128 x 128 x 128 

2.801-11) 

4.00 

2.361-8) 

3.99 

5.251-14) 

5.98 

8 .021-10) 

4.95 

256 x 256 x 256 

1.811-12) 

3.95 

1.501-9) 

3.97 

2.581-13) 

-2.30 

2.531-11) 

4.99 


Example 5 The test Problem 5 can be written as 

d 2 u d 2 u d 2 u . , 

a? + V + a?' /(w> - m0 ' [<UI - 

ii = g(x, y, z), on dQ, 

where f(x, y, z) and g(X y, z ) are determined from the exact solution 

3 3 3 

x 3 yr'z 


(x 2 +y 2 +z 2 ) 


.2x1.5 ’ 


(34) 


which has singularity at the origin and belongs to H 15 E (e is an arbitrary positive constant). 
It follows from the Sobolev embedding theorem that u e C 6-e . 


In the numerical computation, we also use 7 embedded grids with the coarsest grid 
4x4x4, and the corresponding numerical results by the EXCMG method with e = 10 -12 
are listed in table 9 and 10. From table 9 and 10, one can easily find that the results are the 
same as previous examples, i.e., in both L 2 and L°°-norms, the initial guess w/, is a fifth- 
order approximation of the FD solution «*, the FD solution u/, is fourth-order accurate, and 
the numerical gradient Vm/, is also a fourth-order approximation to the exact gradient Vm, 
while the extrapolated solution S/ 7 achieves sixth-order accuracy but starts to lose accuracy 
on the finest grid since a uniform tolerance e = 10 -12 is used on each level of grid. 

We further carry out the computations for other cases when the exact solution u has 
lower regularities, we find that if the exact solution u e H s (s < 7.5), then the extrapolated 
solution u/j will not reach sixth-order accuracy in L°°-norm. In addition, we find that only 
when the exact solution u satisfies that u e H s (s > 5.5), then the numerical solution m/, 
can reach fourth-order accuracy in L“-norm. This is not surprising since H 1 s+e (Q) can 
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Table 11 Errors and convergence rates with e = 10 13 in /.“-norm for Example 6. 


mesh 

IK - “lb 

order 

l|V(«/, - «)||2 

order 

IK - »ll2 

order 

II W h ~ Uh\\l 

order 

16 x 16 x 16 

2.39(—08) 


7.11(—06) 


2.681-09) 


1.541-06) 


32 x 32 x 32 

1.55C-09) 

3.94 

4.45(-07) 

4.00 

6.351-11) 

5.40 

6.71(-08) 

4.52 

64 X 64 X 64 

9.90(—11) 

3.97 

2.46(—08) 

4.18 

1.451-12) 

5.45 

2.391-09) 

4.81 

128 x 128 x 128 

6.26C-12) 

3.98 

1.28C-09) 

4.27 

3.261-14) 

5.48 

7.901-11) 

4.92 

256 x 256 x 256 

3.84(—13) 

4.03 

6.29(-l 1) 

4.34 

1.031-13) 

-1.66 

2.541-12) 

4.96 

Table 12 Errors and 

convergence 

rates with e 

= Kr 13 in L” 

-norm for Example 6. 




mesh 

IK - kIU 

order 

l|V(«/, - w)||oo 

order 

\\Uh w||oo 

order 

II w h - u h \U 

order 

16 x 16 x 16 

1.25(—07) 


3.261-5) 


1.091-07) 


7.051-06) 


32 x 32 x 32 

7.81 (-09) 

4.00 

4.081-6) 

3.00 

6.821-09) 

3.50 

4.401-07) 

4.00 

64 x 64 x 64 

4.88(—10) 

4.00 

5.101-7) 

3.00 

4.261-10) 

3.50 

2.751-08) 

4.00 

128 x 128 x 128 

3.05(—11) 

4.00 

6.381-8) 

3.00 

2.671-11) 

3.50 

1.721-09) 

4.00 

256 x 256 x 256 

2.04(—12) 

3.90 

7.921-9) 

3.01 

3.711-12) 

3.61 

1.081-10) 

4.00 


be continuously embedding into C 6 (Q) and H 5 5+e (Q) can be continuously embedding into 
C 4 (£2) from the Sobolev embedding theorem. 

In the final part of this section, we will show the results for one example where the exact 
solution u e H 5i ~ E (e is an arbitrary small positive constant). 

Example 6 The test Problem 6 can be written as 

d 2 u d 2 u d 2 u Sxyz . _ , 

-T 3-T 3-r = —x-r-— t—— 7 , in Q = [0,1] , 

dx 2 dy 2 dz 2 ( x 2 + y 2 + z 2 ) 0 - 5 (35) 

u = g(x,y,z), ondO, 

where eq. (35) has singularity at the origin and g(x,y, z) is determined from the exact solution 

u = xyzix 1 + y 2 + z 2 ) 0 ' 5 , 

which belongs to H 55 ~ E (e is an arbitrary small positive constant). It follows from the 
Sobolev embedding theorem that u e C 4-e . 

Once again, we use 7 embedded grids with the coarsest grid 4x4x4, and the correspond¬ 
ing numerical results by the EXCMG method with e = 10 -13 are listed in table 11 and 12. 
Since in this case, the exact solution it is only has a finite regularity H 5 5 ~ E (e is any positive 
constant). From table 11 and 12, we can see that the numerical solution u /, is a fourth-order 

approximation to the exact solution in both L 2 and L“-norms. However, due to the lack of 

regularity of the exact solution, we can see from table 11 and 12 that the numerical gradient 
Vm/, converges with fourth-order accuracy in L 2 -norm but only third-order in L°°-norm, the 
extrapolated solution 5/, is 5.5th-order accurate in L 2 -norm but only 3.5th-order accurate in 
L“-norm, while the initial guess w/, is only a fourth-order approximation to the FD solution 
m/, in L°°-norm. Nonetheless, the initial guess w/, is still a fifth-order approximation to the 
FD solution »/, in L 2 -norm. Since the relative residual in the CG solver in our new EXCMG 
method is calculated based on the L 2 -norm (see line 7 of the algorithm 1), thus, our EXCMG 
method is still effective for such low regularity problems (u e H 5 5 ~ E ), and extrapolation can 
also help us to increase the accuracy of initial guess w/, in L 2 -norm, which would widen the 
scope of applicability of our method. 
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Table 13 Comparison of the number of iterations, CPU times (in seconds) and errors between the EXCMG 
method and classical multigrid methods with the Gauss-Seidel smoother. Here CPUw,, denotes the computa¬ 
tional time for constructing of the initial guess w;,. 



e 


V(l,l) 


W(2,l) 



EXCMG 


Iters 1 

CPU 

II u h - u||„ 

Iters 2 

CPU 

\\Ufj m||oo 

Iters J 

CPU 

IK - «IU 

CPU Wft 

Exam 1 

to -' 4 

16 

46.1 

8.61(—13) 

12 

47.6 

8.39(—13) 

8 

12.9 

8.49(-13) 

0.6 

Exam 2 

to - 14 

16 

46.4 

1.91(—12) 

12 

47.6 

1.71 (-12) 

9 

12.6 

1.81(—12) 

0.6 

Exam 3 

10 -“ 

13 

41.5 

7.83(—11) 

9 

39.4 

7.22(— 11) 

8 

11.8 

2.40(-l 1) 

0.6 

Exam 4 

10-° 9 

72 

204.3 

8.27(—10) 

47 

182.9 

2.76(-10) 

8 

10.8 

1.82(—10) 

0.6 

Exam 5 

to - 12 

14 

42.8 

1,77(—12) 

10 

41.3 

1.75(—12) 

9 

13.8 

1.81(—12) 

0.6 

Exam 6 

IQ - 13 

15 

45.9 

1.91(—12) 

11 

46.1 

1.9K-12) 

9 

13.3 

2.04(-12) 

0.6 


1 Iters denotes the number of V-cycles required to reach the error tolerance e. 

2 Iters denotes the number of W-cycles required to reach the error tolerance e. 

3 Iters denotes the number of CG iterations on the finest grid for EXCMG method. 


Table 14 Comparison of the number of iterations, CPU times (in seconds) and errors between the EXCMG 
method and classical multigrid methods with the CG smoother. 



e 


V(l,l) 


W(2,l) 


EXCMG 

Iters 

CPU 

IK - «l|oo 

Iters 

CPU 

IK - «lloo 

Iters 

CPU 

IK - «IU 

Exam 1 

IF" 

15 

43.9 

8.38(—13) 

13 

49.7 

8.42(—13) 

8 

12.9 

8.49(—13) 

Exam 2 

10 14 

15 

43.3 

1.76(—12) 

13 

49.0 

1.73(—12) 

9 

12.6 

1.81(—12) 

Exam 3 

10 -“ 

11 

32.2 

7.09(-ll) 

10 

38.0 

7.22(-l 1) 

9 

11.8 

2.40(-l 1) 

Exam 4 

O 

i 

o 

101 

295.3 

1.74(—10) 

30 

115.8 

1.71(-10) 

8 

10.8 

1.82(—10) 

Exam 5 

10- 12 

13 

39.1 

1.75(-12) 

11 

42.7 

1.75(—12) 

9 

13.8 

1.81(—12) 

Exam 6 

10- 13 

14 

40.9 

1,97(—12) 

11 

42.4 

1.95(—12) 

9 

13.3 

2.04(—12) 


5.2 Computational efficiency 

In this subsection, we compare the efficiency of the our new EXCMG method with the ef¬ 
ficiency of the classical V-cycle and W-cycle multigrid methods for above six examples. 
Results with Gauss-Seidel smoother are listed in table 13 while results with CG smoother 
are listed in table 14. In both tables, the number of iterations, computational time, the L°°- 
norm of the difference between the FD solution it/, and the exact solution u are provided. 
Moreover, the computational time for constructing of the initial guess w/, (line 6 in algo¬ 
rithm 1) is also listed in the last column of table 13, this step contains the extrapolation 
and quartic interpolation as described in section 4.3. By comparing the total computational 
time of the new EXCMG method with the classical V-cycle and W-cycle multigrid methods 
for all above six examples as listed in both table 13 and 14, one can easily find that the 
new EXCMG method needs the smallest time for all examples, and this is particularly true 
when using the unequal meshsizes in different directions, see example 4. Thus, the EXCMG 
method is much more efficient than the classical V-cycle and W-cycle multigrid methods. 
And from the last column in table 13, one can find that the computational time for construct¬ 
ing the initial guess w/, described in section 4.3 is 0.6 seconds for every example, which is 
very cheap, comparing to the total computational time. 

Moreover, one can see from table 13 and 14 that there is only a few number of iterations 
are needed on the finest grid for every example in our EXCMG method, because that the 
initial guess w/, is already an extremely accurate approximation to the FD solution Uh- For 
example, from the last column of table 2, we see that the maximum error of the initial guess 
on the finest grid for example 1 is 4.15 x 10 -12 , which implies that the number of significant 
figures of the approximation exceeds 10. Indeed, from table 1-12 we see that the extrapolated 















An EXCMG method combined with a fourth-order compact scheme for 3D Poisson equation 


19 



4 6 

Iteration 


Fig. 4 Relative residual vs. the number of iterations on the finest grid. 


value Wh served as an initial guess of the CG solver is a fifth-order approximation to the FD 
solution Uh in L 2 -norm, which is one order higher than the convergence order of the fourth- 
order difference solution m/,. Thus, the relative effect of how w/, approximates m;, becomes 
better when mesh is refined, and the number of iterations is reduced most significantly on 
the finest grid, see a more detailed discussion in [44]. 

Finally, we present the curve of the relative residual on the finest grid versus the number 
of iterations for the above six examples in Fig. 4. As we can see that the initial relative 
residual on the finest grid for each example is very small. And due to the high oscillations 
of the initial error as shown in section 4.4, the relative residual decreases by several orders 
of magnitude after only a few iterations, and then reaches a number that is less than the 
required tolerance. 


6 Conclusions 

In this work, we developed a new extrapolation cascadic multigrid (EXCMG) method com¬ 
bined with 19-point fourth-order compact difference scheme for solving the 3D Poisson 
equation on rectangular domains. The major advantage of the method is to use the Richard¬ 
son extrapolation and tri-quartic Lagrange interpolation techniques for two numerical solu¬ 
tions on two-level of grids (current and previous grids) to obtain a fifth-order approximation 
Wh to the fourth-order FD solution Uh as the initial guess of the iterative solution on the next 
finer grid, which greatly reduces the iteration numbers. When the exact solution u is suffi¬ 
ciently smooth, a sixth-order extrapolated solution S/, on the fine grid can be obtained by 
using two fourth-order numerical solutions on two scale grids. Moreover, the gradient of so¬ 
lution Wit/, can also be computed easily and efficiently through solving a series of tridiagonal 
linear systems resulting from the fourth-order compact FD discretization of the derivatives. 
Finally, numerical results show that our new extrapolation cascadic multigrid method is 
much more efficient comparing to the classical V-cycle and W-cycle multigrid method and 
it is particularly suitable for solving large scale problems. 

The work presented in this paper is an extension of our previous work, which is based 
on the EXCMG method for the 3D elliptic problem with the linear FE discretization [44]. 
In the near future, we will extend our method to convection-diffusion equations, Flelmholtz 
equations, biharmonic equations, and other related equations. 























20 


Kejia Pan et al. 


Acknowledgements Kejia Pan was supported by the National Natural Science Foundation of China (Nos. 
41474103 and 41204082), the National High Technology Research and Development Program of China (No. 
2014AA06A602), the Natural Science Foundation of Hunan Province of China (No. 2015JJ3148). Dongdong 
He was supported by the National Natural Science Foundation of China (No. 11402174), the Program for 
Young Excellent Talents at Tongji University (No. 2013KJ012) and the Scientific Research Foundation for the 
Returned Overseas Chinese Scholars, State Education Ministry. Hongling Hu was supported by the National 
Natural Science Foundation of China (No. 11301176). 


References 

1. J.C. Strikwerda, Finite difference schemes and partial differential equations. Chapman & Hall, 1989. 

2. M.M. Gupta, A fourth-order Poisson solver, J. Comput. Phys. 55 (1985) 166-172. 

3. M.M. Gupta, J. Kouatchou, Symbolic derivation of finite difference approximations for the three- 
dimensional Poisson equation, Numer. Meth. Part. D. E. 14 (1998) 593-606. 

4. W.F. Spotz, G.F. Carey, A high-order compact formulation for the 3D poisson equation, Numer. Meth. 
Part. D. E. 12 (1996) 235-243. 

5. G. Sutmann, B. Steffen, High-order compact solvers for the three-dimensional Poisson equation, J. Com¬ 
put. Appl. Math. 187 (2006) 142-170. 

6 . J. Wang, W. Zhong, J. Zhang, A general meshsize fourth-order compact difference discretization scheme 
for 3D Poisson equation, Appl. Math. Comput. 183 (2006) 804-812. 

7. M.M. Gupta, J. Kouatchou, J. Zhang, Comparison of second-order and fourth-order discretization for 
multigrid Poisson solvers, J. Comput. Phys. 132 (1997) 226-232. 

8 . M. Othman, A.R. Abdullah, An efficient multigrid Poisson solver, Int. J. Comput. Math. 71 (1999) 541- 
553. 

9. S. Schaffer, High order multi-grid methods, Math. Comp. 43 (1984) 89-115. 

10. J. Zhang, Multigrid method and fourth-order compact scheme for 2D Poisson equation with unequal 
mesh-size discretization, J. Comput. Phys. 179 (2002) 170-179. 

11. Y. Wang, J. Zhang, Sixth-order compact scheme combined with multigrid method and extrapolation 
technique for 2D poisson equation, J. Comput. Phys. 228 (2009) 137-146. 

12. J. Zhang, Fast and high accuracy multigrid solution of the three dimensional Poisson equation, J. Comput. 
Phys. 143 (1998) 449-161. 

13. Y.B. Ge, Multigrid method and fourth-order compact difference discretization scheme with unequal 
meshsizes for 3D poisson equation, J. Comput. Phys. 229 (2010) 6381-6391. 

14. S.F. McCormick (Ed.), Multigrid Methods, Frontiers in Applied Mathematics, SIAM, Philadelphia, PA, 
1987. 

15. W.L. Briggs, S.F. McCormick, V. E. Henson, A Multigrid Tutorial, second ed., SIAM, Philadelphia, PA, 
2000. 

16. U. Trottenberg, C.W. Oosterlee, A. Schller, Multigrid, Academic Press, London, 2001. 

17. H. Moghaderi, M. Dehghan, M. Hajarian, A fast and efficient two-grid method for solving d-dimensional 
Poisson equations, Numer. Algor. 72 (2016) 483C537. 

18. I. Altas, J. Dym, M.M. Gupta, R.P. Manohar, Multigrid solution of automatically generated high-order 
discretizations for the biharmonic equation, SIAM J. Sci. Comput. 19 (1998) 1575-1585. 

19. J. Zhang, H. Sun, J.J. Zhao, High order compact scheme with multigrid local mesh refinement procedure 
for convection diffusion problems, Comput. Method. Appl. M. 191 (2002) 4661-4674. 

20. Y.B. Ge, F.J. Cao, Multigrid method based on the transformation-free HOC scheme on nonuniform grids 
for 2D convection diffusion problems, J. Comput. Phys. 230 (2011) 4051-4070. 

21. Y. Wang, J. Zhang, Fast and robust sixth-order multigrid computation for the three-dimensional 
convection-diffusion equation, J. Comput. Appl. Math. 234 (2010) 3496-3506. 

22. F.A. Bomemann, P. Deuflhard, The cascadic multigrid method for elliptic problems, Numer. Math. 75 
(1996) 135-152. 

23. V. Shaidurov, Some estimates of the rate of convergence for the cascadic conjugate-gradient method, 
Comput. Math. Appl. 31 (1996) 161-171. 

24. D. Braess, W. Dahmen, A cascadic multigrid algorithm for the stokes equations, Numer. Math. 82 (1999) 
179-191. 

25. G. Timmermann, A cascadic multigrid algorithm for semilinear elliptic problems, Numer. Math. 86 
(2000)717-731. 

26. V. Shaidurov, L. Tobiska, The convergence of the cascadic conjugate- gradient method applied to elliptic 
problems in domains with re-entrant cor- ners, Math. Comput. 69 (2000) 501-520. 



An EXCMG method combined with a fourth-order compact scheme for 3D Poisson equation 


21 


27. V. Shaidurov, G. Timmermann, A cascadic multigrid algorithm for semi- linear indefinite elliptic prob¬ 
lems, Computing 64 (2000) 349-366. 

28. Z.C. Shi, X.J. Xu, Cascadic multigrid for parabolic problems, J. Comput. Math. 18 (2000) 551-560. 

29. D. Braess, P. Deuflhard, K. Lipnikov, A subspace cascadic multigrid method for mortar elements, Com¬ 
puting 69 (2002) 205-225. 

30. R. Stevenson, Nonconforming finite elements and the cascadic multi-grid method, Numer. Math. 91 
(2002) 351-387. 

31. S.Z. Zhou, H.X. Hu, On the convergence of a cascadic multigrid method for semilinear elliptic problem, 
Appl. Math. Comput. 159 (2004) 407417. 

32. Q. Du, P.B. Ming, Cascadic multigrid methods for parabolic problems, Sci. China Ser. A-Math. 51 (2008) 
1415-1439. 

33. X.J. Xu, W.B. Chen, Standard and economical cascadic multigrid methods for the mortar finite element 
methods, Numer. Math-Theory Me. 2 (2009) 180-201. 

34. H.X. Yu, J.P. Zeng, A cascadic multigrid method for a kind of semilinear elliptic problem, Numer. Al- 
gorightm 58 (2011) 143-162. 

35. Z.C. Shi, X.J. Xu, Y.Q. Huang. Economical cascadic multigrid method (ECMG). Science in China Series 
A: Mathematics. 2007, 50(12): 1765 -1780. 

36. C.M. Chen, H.L. Hu, Z.Q. Xie, et al, Analysis of extrapolation cascadic multigrid method (EXCMG), 
Sci. China Ser. A-Math. 51 (2008) 1349-1360. 

37. C.M. Chen, Z.C. Shi, H.L. Hu, On extrapolation cascadic multigrid method, J. Comput. Math. 29 (2011) 
684-697. 

38. H.L. Hu, C.M. Chen, K.J. Pan, Asymptotic expansions of finite element solutions to Robin problems in 
H 3 and their application in extrapolation cascadic multigrid method, Sci. China Math. 57 (2014) 687- 
698. 

39. H.L. Hu, C.M. Chen, K.J. Pan, Time-extrapolation algorithm (TEA) for linear parabolic problems, J. 
Comput. Math. 32 (2014) 183-194. 

40. K.J. Pan, J.T. Tang, H.L. Hu, et al, Extrapolation cascadic multigrid method for 2.5D direct current 
resistivity modeling (in Chinese), Chinese J. Geophys. 55 (2012) 2769-2778. 

41. K.J. Pan, J.T. Tang, 2.5-D and 3-D DC resistivity modelling using an extrapolation cascadic multigrid 
method, Geophys. J. Int. 197 (2014) 1459-1470. 

42. G.A. Newman, A Review of high-performance computational strategies for modeling and imaging of 
electromagnetic induction data, Surv. Geophys. 35 (2014) 85-100. 

43. G. Berikelashuili, M.M. Gupta, M. Mirianashvili, Convergence of fourth-order compact difference 
schemes for three-dimensional convection-diffusion equations, SIAM J. Numer. Anal. 45 (2007) 443- 
455. 

44. K.J. Pan, D.D. He, H.L. Hu, A new extrapolation cascadic multigrid method for 3D elliptic boundary 
value problems on rectangular domains, arXiv preprint arXiv: 1506.02983 (2015). 

45. G.I. Marchuk, V.V. Shaidurov, Difference Methods and Their Extrapolations, Springer-Verlag, New 
York, 1983. 

46. P. Neittaanmaki, Q. Lin, Acceleration of the convergence in finite-difference method by predictor cor¬ 
rector and splitting extrapolation methods, J. Comput. Math. 5 (1987) 181-190. 

47. FoBmeier R, On Richardson extrapolation for finite difference methods on regular grids, Num. Math. 55 
(1989)451-462. 

48. G.Q. Han, Spline finite difference methods and their extrapolation for singular two-point boundary value 
problems, J. Comput. Math. 11 (1993) 289-296. 

49. H. Sun, J. Zhang, A high order finite difference discretization strategy based on extrapolation for con¬ 
vection diffusion equations, Numer. Meth. Part. D. E. 20 (2004) 18-32. 

50. K. Rahul, S.N. Bhattacharyya, One-sided finite-difference approximations suitable for use with Richard¬ 
son extrapolation, J. Comput. Phys. 219 (2006) 13-20. 

51. J.B. Munyakazi, K.C. Patidar, On Richardson extrapolation for fitted operator finite difference methods, 
Appl. Math. Comput. 201 (2008) 465-480. 

52. C.K.W. Tam, K.A. Kurbatskii, A wavenumber based extrapolation and interpolation method for use in 
conjunction with high-order finite difference schemes, J. Comput. Phys. 157 (2000) 588-617. 

53. Y. Ma, Y. Ge, A high order finite difference method with Richardson extrapolation for 3D convection 
diffusion equation, Appl. Math. Comput. 215 (2010) 3408-3417. 

54. C.H. Marchi, L.A. Novak, C.D. Santiago, et al, Highly accurate numerical solutions with repeated 
Richardson extrapolation for 2D Laplace equation, Appl. Math. Model. 37 (2013) 7386-7397. 

55. L. Collatz. The Numerical Treatment of Differential Equations, Springer-Verlag, New York, 1966, P.538. 



